Dataload and process

First part

data filter and normalization

source("./utils/tianfengRwrappers.R")

Loading required package: reticulate
Loading required package: tidyr

Attaching package: ‘tidyr’

The following object is masked from ‘package:S4Vectors’:

    expand


Attaching package: ‘MySeuratWrappers’

The following objects are masked from ‘package:Seurat’:

    DimPlot, DoHeatmap, LabelClusters, RidgePlot, VlnPlot


Attaching package: ‘cowplot’

The following object is masked from ‘package:ggpubr’:

    get_legend

Loading required package: viridisLite

Attaching package: ‘reshape2’

The following object is masked from ‘package:tidyr’:

    smiths

NOTE: Either Arial Narrow or Roboto Condensed fonts are required to use these themes.
      Please use hrbrthemes::import_roboto_condensed() to install Roboto Condensed and
      if Arial Narrow is not on your system, please see https://bit.ly/arialnarrow

clusterProfiler v3.16.1  For help: https://guangchuangyu.github.io/software/clusterProfiler

If you use clusterProfiler in published research, please cite:
Guangchuang Yu, Li-Gen Wang, Yanyan Han, Qing-Yu He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology. 2012, 16(5):284-287.

Attaching package: ‘clusterProfiler’

The following object is masked from ‘package:DelayedArray’:

    simplify

The following object is masked from ‘package:IRanges’:

    slice

The following object is masked from ‘package:S4Vectors’:

    rename

The following object is masked from ‘package:stats’:

    filter

Registering fonts with R

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:IRanges’:

    slice

The following object is masked from ‘package:S4Vectors’:

    rename

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

Loading required package: e1071
Warning: couldn't connect to display ":0"
Attaching package: ‘widgetTools’

The following object is masked from ‘package:dplyr’:

    funs


Attaching package: ‘DynDoc’

The following object is masked from ‘package:DelayedArray’:

    path

The following object is masked from ‘package:BiocGenerics’:

    path


Attaching package: ‘IDPmisc’

The following object is masked from ‘package:shiny’:

    hr


Attaching package: ‘DT’

The following objects are masked from ‘package:shiny’:

    dataTableOutput, renderDataTable

The following object is masked from ‘package:SeuratObject’:

    JS

The following object is masked from ‘package:Seurat’:

    JS

========================================
circlize version 0.4.15
CRAN page: https://cran.r-project.org/package=circlize
Github page: https://github.com/jokergoo/circlize
Documentation: https://jokergoo.github.io/circlize_book/book/

If you use it in published research, please cite:
Gu, Z. circlize implements and enhances circular visualization
  in R. Bioinformatics 2014.

This message can be suppressed by:
  suppressPackageStartupMessages(library(circlize))
========================================

Loading required package: grid
========================================
ComplexHeatmap version 2.4.3
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference

If you use it in published research, please cite:
Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
  genomic data. Bioinformatics 2016.

This message can be suppressed by:
  suppressPackageStartupMessages(library(ComplexHeatmap))
========================================


Attaching package: ‘ComplexHeatmap’

The following object is masked from ‘package:plotly’:

    add_heatmap

Loading required package: AnnotationDbi

Attaching package: ‘AnnotationDbi’

The following object is masked from ‘package:plotly’:

    select

The following object is masked from ‘package:clusterProfiler’:

    select

The following object is masked from ‘package:dplyr’:

    select

From cellranger

QC

fly_merge <- CreateSeuratObject(Read10X("./raw_data_file/fly_merge/"), names.field = 2, names.delim = "-",
                                   project = "fly_merge", min.cells = 10, min.features = 200)

fly_merge <- PercentageFeatureSet(fly_merge, pattern = "^mt:", col.name = "percent.mt")

VlnPlot(fly_merge,c("percent.mt"))/
VlnPlot(fly_merge,c("nFeature_RNA"))/
VlnPlot(fly_merge,c("nCount_RNA"))

fly_merge <- subset(fly_merge, subset = nFeature_RNA > 300 & nFeature_RNA < 4000 & 
        nCount_RNA > 1000 &  nCount_RNA < 20000 & percent.mt < 10)

fly_merge <- fly_merge %>% SCTransform(vars.to.regress = "percent.mt", verbose = F) %>% 
  RunPCA() %>% FindNeighbors(dims = 1:20) %>% 
  RunUMAP(dims = 1:20) %>% 
  FindClusters(resolution = 0.1)

PC_ 1 
Positive:  regucalcin, CG31777, CG8501, CG5399, Gal, CG3940, CG14629, Gs1, Glt, CG4250 
       fat-spondin, Arpc3B, cer, lectin-28C, Hsp23, Ama, CG44250, Npc2a, CG31337, CG42566 
       Col4a1, scf, CG34296, fax, CG4950, Hexo2, SPARC, firl, CG10131, lectin-24Db 
Negative:  Phk-3, Gasp, ImpE2, ImpE1, Spn43Aa, CG15353, ImpL1, hui, CG2852, Spn100A 
       CG12009, Antp, obst-B, eEF1alpha1, RpL41, CG2016, serp, Fas2, CG9691, hth 
       CG2663, dsx-c73A, Fas3, pio, ogre, CG31997, CG11370, Cpr51A, Cyp1, betaTub56D 
PC_ 2 
Positive:  Gasp, Spn43Aa, CG12009, Cpr51A, Osi15, CG15353, CG11370, obst-B, CG7220, CG2016 
       CAH2, dpy, serp, l(2)34Fc, mt:lrRNA, Tsf1, peb, CG30154, sn, betaTub60D 
       CG8854, uif, Osi20, CG30197, form3, CG9782, obst-A, vvl, CG2852, lncRNA:CR44811 
Negative:  Phk-3, ImpL1, ImpE1, ImpE2, dsx-c73A, Antp, CG13737, Spn100A, CG34232, GstD1 
       hui, CG2663, CG16798, ct, CG32354, vir-1, CG9691, qtc, Ect3, Cpr47Eb 
       frm, Mur89F, Oaz, Cad88C, CG31637, tyn, slbo, lncRNA:CR43432, Cyp18a1, tnc 
PC_ 3 
Positive:  Gasp, CG2663, Spn100A, Cpr51A, Osi15, CG34232, CG13737, CG4702, CG16798, CG12009 
       dsx-c73A, ect, CG15353, Phk-3, qtc, CG32354, Ect3, CG2016, vir-1, Cad88C 
       Spn43Aa, CG7220, Antp, lncRNA:CR43314, Nep2, Cpr47Eb, Ldh, obst-A, CAH2, CG2852 
Negative:  Him, E(spl)m6-BFM, E(spl)m7-HLH, E(spl)mbeta-HLH, CG11835, E(spl)malpha-BFM, E(spl)m2-BFM, twi, E(spl)m8-HLH, hdc 
       hth, eEF1alpha2, E(spl)mgamma-HLH, ImpE2, vg, E(spl)m3-HLH, CG9650, Ppn, stg, E(spl)mdelta-HLH 
       side, Act57B, zfh1, Act87E, Pen, beat-IIIc, bib, smal, Hsp26, Sp1 
PC_ 4 
Positive:  CG34232, CG2663, CG13737, Him, E(spl)m6-BFM, qtc, ect, CG16798, Cad88C, CG32354 
       Osi15, Ect3, Cpr51A, CG4702, Gasp, dsx-c73A, slbo, Antp, CG10737, E(spl)m7-HLH 
       CG11835, Cpr47Eb, bol, SPARC, cpo, Oaz, Ldh, CG7724, twi, betaTub97EF 
Negative:  Phk-3, CG15353, Idgf4, CG9691, GstD1, ImpE2, Mmp1, Nep2, ImpE1, CG11370 
       Reg-5, l(2)34Fc, tnc, CG4928, ple, Cpr12A, frm, Spn43Aa, b, AANAT1 
       Cpr49Ac, ImpL1, CG9312, CG9192, Tsf1, Dl, Ance, Gs2, CG6770, Ddc 
PC_ 5 
Positive:  CG15353, CG11370, Spn100A, CG2663, CG34232, l(2)34Fc, Tsf1, Nep2, Spn43Aa, CG2016 
       CG13737, Gs2, lncRNA:CR43314, Ddc, l(2)k05911, Idgf4, qtc, CG13748, CG7888, Mmp1 
       CG13639, CG13082, bnl, obst-B, CG16798, CG3502, CG32354, Cad88C, SP1029, obst-A 
Negative:  Phk-3, ImpE2, Gasp, Cpr51A, ImpL1, Osi15, ImpE1, GstD1, CG7220, CG12009 
       Osi20, hui, Dbi, lncRNA:CR43432, lncRNA:CR44811, betaTub60D, CG9689, pio, bs, CG4702 
       CG17121, CG1368, CG15887, tsr, mirr, sn, frm, Mur89F, Rtnl1, chic 
Computing nearest neighbor graph

Computing SNN
Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session10:33:01 UMAP embedding parameters a = 0.9922 b = 1.112
10:33:01 Read 18195 rows and found 20 numeric columns
10:33:01 Using Annoy for neighbor search, n_neighbors = 30
10:33:01 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
10:33:04 Writing NN index file to temp file /tmp/RtmpLc7AED/file1099d32841635
10:33:04 Searching Annoy index using 1 thread, search_k = 3000
10:33:09 Annoy recall = 100%
10:33:10 Commencing smooth kNN distance calibration using 1 thread
10:33:12 Initializing from normalized Laplacian + noise
10:33:13 Commencing optimization for 200 epochs, with 769526 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
10:33:21 Optimization finished

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 18195
Number of edges: 596302

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9703
Number of communities: 11
Elapsed time: 1 seconds

table(fly_merge$orig.ident)

    1     2     3 
10191   989  7015 

trachea <- subset(fly_merge, idents = c(0:3,6:8))
trachea <- trachea %>% 
  RunPCA() %>% FindNeighbors(dims = 1:20) %>% 
  RunUMAP(dims = 1:20)
DimPlot(trachea, cells.highlight = WhichCells(larva))
trachea <- trachea %>% FindClusters(resolution = 0.2)

trachea$ctrl_idents <- larva$seurat_clusters
fly_merge$ctrl_idents <- larva$seurat_clusters

ctrl <- subset(fly_merge_filtered, orig.ident == 1)
# saveRDS(ctrl,"./processed_data_file/ctrl.rds")

HSD <- subset(fly_merge_filtered, orig.ident == 2 | orig.ident == 3)
# saveRDS(HSD,"./processed_data_file/HSD.rds")

run harmony

library(harmony)
umapplot(trachea, split.by = "orig.ident")


trachea_harmony <- RunHarmony(trachea, "orig.ident", plot_convergence = T, project.dim = F, assay.use = "SCT")
trachea_harmony@reductions[["pca"]] <- trachea_harmony@reductions[["harmony"]]


trachea_harmony <- trachea_harmony %>% RunUMAP(dims = 1:20)
trachea_harmony <- FindClusters(trachea_harmony, resolution = 0.25)
umapplot(trachea_harmony)
umapplot(trachea_harmony, group.by = "ctrl_idents")
# saveRDS(fly_merge_harmony,"./processed_data_file/fly_merge_harmony.rds") #已经经过分组处理了

ctrl_filtered <- NormalizeData(ctrl_filtered, assay = "RNA") %>% ScaleData(assay = "RNA")

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Centering and scaling data matrix

  |                                                                                       
  |                                                                                 |   0%
  |                                                                                       
  |=======                                                                          |   8%
  |                                                                                       
  |==============                                                                   |  17%
  |                                                                                       
  |====================                                                             |  25%
  |                                                                                       
  |===========================                                                      |  33%
  |                                                                                       
  |==================================                                               |  42%
  |                                                                                       
  |========================================                                         |  50%
  |                                                                                       
  |===============================================                                  |  58%
  |                                                                                       
  |======================================================                           |  67%
  |                                                                                       
  |=============================================================                    |  75%
  |                                                                                       
  |====================================================================             |  83%
  |                                                                                       
  |==========================================================================       |  92%
  |                                                                                       
  |=================================================================================| 100%

HSD_filtered <- NormalizeData(HSD_filtered, assay = "RNA") %>% ScaleData(assay = "RNA")

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Centering and scaling data matrix

  |                                                                                       
  |                                                                                 |   0%
  |                                                                                       
  |=======                                                                          |   8%
  |                                                                                       
  |==============                                                                   |  17%
  |                                                                                       
  |====================                                                             |  25%
  |                                                                                       
  |===========================                                                      |  33%
  |                                                                                       
  |==================================                                               |  42%
  |                                                                                       
  |========================================                                         |  50%
  |                                                                                       
  |===============================================                                  |  58%
  |                                                                                       
  |======================================================                           |  67%
  |                                                                                       
  |=============================================================                    |  75%
  |                                                                                       
  |====================================================================             |  83%
  |                                                                                       
  |==========================================================================       |  92%
  |                                                                                       
  |=================================================================================| 100%

END

celltype annotation